{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install dependencies for this example\n",
    "# Note: This does not include itkwidgets, itself\n",
    "import sys\n",
    "!{sys.executable} -m pip install --upgrade itk-minimalpathextraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import itk\n",
    "import numpy as np\n",
    "\n",
    "from itkwidgets import view"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_DIR = Path('../test/Input')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0072c3bbb27a4bd4b35a965458611254",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[], gradient_opacity=0.22, point_sets=[], rendered_image=<itkImagePython.itkImageUC2; proxy …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dsa = itk.imread(str(DATA_DIR / 'Real-DSA-01.jpg'))\n",
    "view(dsa, ui_collapsed=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6e9289c05fec4f19b77ebab14a5267e4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[], gradient_opacity=0.22, point_sets=[], rendered_image=<itkImagePython.itkImageF2; proxy o…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dsa_speed = itk.imread(str(DATA_DIR / 'Real-DSA-01-Speed-02.mhd'))\n",
    "view(dsa_speed, ui_collapsed=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Path: [394.00, 203.00] [388.00, 229.00] [356.00, 217.00] [289.00, 252.00] [194.00, 309.00] [162.00, 318.00] [232.00, 490.00]\r\n"
     ]
    }
   ],
   "source": [
    "path_file = str(DATA_DIR / 'Real-DSA-01.path')\n",
    "!cat {path_file}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "PointType = itk.Point[itk.D,2]\n",
    "PathInformationType = itk.SpeedFunctionPathInformation[PointType]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_information = PathInformationType.New()\n",
    "with open(path_file, 'r') as fp:\n",
    "    for line in fp:\n",
    "        line = line.replace(\"Path: \", \"\")\n",
    "        line = line.replace(\"[\", \"\").strip()\n",
    "        points = line.split(']')[:-1]\n",
    "        way_points = np.empty((0, 2))\n",
    "        for index, point in enumerate(points):\n",
    "            point_float = PointType()\n",
    "            point_float[0] = float(point.split(',')[0])\n",
    "            point_float[1] = float(point.split(',')[1])\n",
    "            if index == 0:\n",
    "                start_point = np.array(point_float).reshape(1,2)\n",
    "                path_information.SetStartPoint(point_float)\n",
    "            elif index == len(points)-1:\n",
    "                end_point = np.array(point_float).reshape(1,2)\n",
    "                path_information.SetEndPoint(point_float)\n",
    "            else:\n",
    "                way_points = np.vstack((way_points, np.array(point_float).reshape(1,2)))\n",
    "                path_information.AddWayPoint(point_float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d907dae42c9243cf862bbb75e3f2ba77",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[], gradient_opacity=0.22, point_set_colors=array([[1. , 0. , 0. ],\n",
       "       [1. , 0. , 0.5],\n",
       "…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view(dsa_speed,\n",
    "     point_sets=[start_point, way_points, end_point],\n",
    "     point_set_colors=[(1.0,0.0,0.0), (1.0,0.0,0.5), (1.0,0.0,0.8)],\n",
    "     ui_collapsed=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "interpolator = itk.LinearInterpolateImageFunction.New(dsa_speed)\n",
    "interpolator.SetInputImage(dsa_speed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "cost_function = itk.SingleImageCostFunction[type(dsa_speed)].New(interpolator=interpolator)\n",
    "cost_function.SetInterpolator(interpolator)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "spacing = list(dsa_speed.GetSpacing())\n",
    "min_spacing = min(spacing)\n",
    "step_length_factor = 0.25\n",
    "step_length_relax = 0.5\n",
    "number_of_iterations = 4000\n",
    "\n",
    "optimizer = itk.RegularStepGradientDescentOptimizer.New(\n",
    "    number_of_iterations=number_of_iterations,\n",
    "    maximum_step_length=1.0*step_length_factor*min_spacing,\n",
    "    minimum_step_length=0.5*step_length_factor*min_spacing,\n",
    "    relaxation_factor=step_length_relax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "termination_value = 3.0\n",
    "path_filter = itk.SpeedFunctionToPathFilter.New(dsa_speed,\n",
    "                                               cost_function=cost_function,\n",
    "                                               optimizer=optimizer,\n",
    "                                               termination_value=termination_value)\n",
    "path_filter.SetInput(dsa_speed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_filter.AddPathInformation(path_information)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 64.8 ms, sys: 0 ns, total: 64.8 ms\n",
      "Wall time: 64.2 ms\n"
     ]
    }
   ],
   "source": [
    "%time path_filter.Update()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2232"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "path = path_filter.GetOutput(0)\n",
    "path.GetVertexList().Size()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6ed292f79b4647408ba500afa897af1a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'name': '_points', 'numberO…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view(dsa, geometries=path, ui_collapsed=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
